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It is sometimes speculated that the sign problem that afflicts many quantum field theories 
might be reduced or even eliminated by choosing an alternative domain of integration within 
a complexified extension of the path integral (in the spirit of the stationary phase integration 
method). In this paper we start to explore this possibility somewhat systematically. A 
first inspection reveals the presence of many difficulties but — quite surprisingly — most of 
them have an interesting solution. In particular, it is possible to regularize the lattice 
theory on a Lefschetz thimble, where the imaginary part of the action is constant and 
disappears from all observables. This regularization can be justified in terms of symmetries 
and perturbation theory. Moreover, it is possible to design a Monte Carlo algorithm that 
samples the configurations in the thimble. This is done by simulating, effectively, a five 
dimensional system. We describe the algorithm in detail and analyze its expected cost and 
stability. Unfortunately, the measure term also produces a phase which is not constant and 
it is currently very expensive to compute. This residual sign problem is expected to be much 
milder, as the dominant part of the integral is not affected, but we have still no convincing 
evidence of this. However, the main goal of this paper is to introduce a new approach to 
the sign problem, that seems to offer much room for improvements. An appealing feature of 
this approach is its generality. It is illustrated first in the simple case of a scalar field theory 
with chemical potential, and then extended to the more challenging case of QCD at finite 
baryonic density. 



INTRODUCTION 



Formidable experimental efforts are presently being devoted to study strongly interacting nu- 
clear matter in a high density medium (see jl[ for a recent and comprehensive review). In fact, 
quantum chromodynamics (QCD) is expected to display a very rich phase structure in that regime 
But, unfortunately, lattice QCD calculations are severely limited by the lack of a positive 
measure, which is necessary for the direct applicability of importance sampling Monte Carlo meth- 
ods. As a consequence of this sign problem, one expects 0] the cost of direct Monte Carlo methods 
to scale like e v (where V is the four dimensional volume), which is clearly prohibitive. 

In the past decade, much progress has been achieved in devising techniques to alleviate the sign 
problem. The present status has been reviewed in the recent lattice conferences 0-0] ■ Thanks 
to the reweightingmethod jf|, the Taylor expansion method 0, S] and the imaginary chemical 
potential method [9|, [10] it is now possible to perform quite reliable calculations in the important 
region near the finite temperature phase transition at small chemical potential. 

In the high density region, other approaches based on the complex Langevin equation 
13], on worm algorithms on an effective 3d theory 15 . [lrj ]. on the histogram method 
on the factorization or density of state method [18H21I]. on the generalized imaginary chemical 
potential method [22]], on the fugacity expansion [23|], on dimensional reduction 24] and the large 
N c limit 0, We been proposed and present" promising aspects. However, no method has 
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yet demonstrated reliability in the high density regime of QCD. In this context, the search for 
alternatives is certainly very desirable. 

The approach proposed in this paper draws its inspiration from the simple idea of saddle-point 
integration along the paths of steepest descent. This is a powerful and elementary tool to treat 
oscillating, low dimensional, integrals, such as: 



X= [ dxg(x) e /(x) . fl ,/:Hol(C). 

JTSL 



It is important to distinguish two independent conceptual steps in the classic saddle-point integra- 
tion method. The first step consists in deforming the domain of integration R into a path 7 C C 
that preserves the homology class of the original integral. The path 7 typically goes through one 
stationary point of / and then follows the direction of steepest descent of the real part of the action. 
By holomorphicity of /, the imaginary part of / is constant along 7, which justifies the name of 
stationary phase method. The second step consists in Taylor expanding / around the stationary 
point, which typically provides a satisfactory approximation of X. 

This technique is very effective for computing low dimensional oscillating integrals. But, the 
approximation in the second step is not satisfactory for our goal of a non-perturbative formulation 
of a quantum field theory (QFT), as it would amount to some form of perturbative expansion. 
On the other hand, the deformation of the domain of integration, described in the first step is 
potentially very interesting, in relation to the sign problem of QCD. In fact, the imaginary part 
of f(x) is constant along the path 7, and can be factorized out of the integral in the form of a 
constant phase factor. Moreover, the real part of f(x) is bounded from above by its value at the 
stationary point. 

The extension of these ideas to the multi-dimensional case is a classic topic in complex analysis 



271 . |28( |. In this case, the integral to solve has the form: 

1 = / dxi . . . dx n g(x) , 



where the functions / and g are now holomorphic in C n — > C. Under suitable conditions on / and 



g, Picard-Lefschetz/Morse theory [29l . |30| shows that one can associate to each stationary point 
p a G C n (a G S) of the function / an integration domain J a of real dimension n (an n-cycle) 
immersed in C n . Such n-cycles (called Lefschetz thimbles) generalize the idea of path of steepest 
descent and, altogether, they provide a basis of the relevant homology group, so that any n-cycle C, 
along which we might want to integrate g{x)e^ x \ can be expressed in terms of the basis {Jo-jo-es, 



with integer coefficients n a [23], i.e.: 



C = n *J* C 1 ) 
o-es 

In the case of QCD at finite baryonic density (QCD/i), the usual functional integral is well 
defined, on a finite lattice, on the integration domain C = SU(3) Vx4 . The latter is an n-cycle 
(where n = 8 x 4 x V) immersed in SL(3,C) Vx4: , that belongs to a well defined homology class, 
that can be also written in terms of the basis of thimbles {Ja} introduced above, with well defined 
coefficients n a . So, in principle, the original functional integral could be expressed as the sum of 
integrals on the thimbles J a , where the integrand has no sign problem (by a generalization of the 
stationary phase property, as we will see). But, finding the stationary points p a , computing the 
coefficients n a 7^ and performing simulations on each J a is not feasible. However, we should also 
ask whether it is necessary. In fact, although the QCD^ partition function already has a well defined 
regularization, on a finite lattice, it might be worthwhile to consider an alternative one, if the latter 
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had some practical advantage. If we adopt this point of view, then our guiding principle should 
be the necessity of constructing a local QFT that reproduces the correct symmetries (including 
the correct representations and degrees of freedom) and the correct perturbative expansion. By 
universality, we expect that these properties essentially determine the scaling behavior in the 
continuum. If so, it is not necessary that the new integration domain belong exactly to the same 
homology class of the original integral, and it is natural to define a lattice regularization of a 
QFT as a functional integral over that single thimble Jo which is associated to the perturbative 
stationary point. Our first task is to show that this regularization has indeed the correct symmetry 
representations and the correct perturbative expansion. 

Recently, Witten [3l| used Morse theory to extend the definition of the three dimensional Chern- 
Simon QFT to a set of complex values of the parameters where the integral is not manifestly 
convergent. In particular, Witten computed analytically how the (unnormalized) Chern-Simon 
partition function depends on the parameters of the action (in presence of a knot background). In 
order to compute such dependence exactly, it was necessary to determine how the coefficients n a of 
Eq. ([1]) depend on the parameters of the action, and, to do that, the so-called Stokes phenomena 
had to be taken into account [311 ]. This is not realistic in the case of QCD^, neither analytically nor 
numerically. However, Monte Carlo methods suggest a different approach. Although the knowledge 
of the parameter dependence of the partition function is certainly the classic and convenient way 
to compute the corresponding observables analytically, it is not the only way. In particular, this is 
never done in a Monte Carlo calculation, where, instead, one fixes the normalization by computing 
a suitable set of observables. Consequently, a uniform normalization of the partition function for 
different values of the parameters is not necessary: the normalization can (and often must) be 
performed independently in each point of the parameter space. From this point of view, it is more 
convenient and natural to regularize the theory always on that thimble Jo which has the correct 
perturbative limit, as suggested above. 

Once we have defined our regularization and proved its properties, we turn to the question of 
how to simulate numerically what we have proposed. This is far from straightforward. In fact, 
although the Lefschetz thimble Jo is a smooth manifold [3Q], it is not clear how to compute the 
tangent space of Jo at a given gauge configuration A = {A®(x)} £ Jo, by using only information 
available in the neighborhood of A. If we adopt the Langevin algorithm to simulate the system, the 
problem is partially solved, because — as we will see — the force term that appears in the Langevin 
equation is tangent to Jo by construction. But this solves the problem only in part, because the 
Langevin dynamics also requires a noise term, and this needs to be projected on the tangent space 
of Jo. However, the tangent space Ta(Jo) is easy to compute at the stationary point A = and, 
for any other configuration A E Jo, we can use the flow of steepest descent to parallel transport 
a tangent vector r\ £ To{Jo) into a vector rf £ Ta(Jo)- The concept of Lie derivative represents 
the natural tool to accomplish the parallel transport along the flow of steepest descent. This 
method presents challenging aspects, from the numerical point of view, but we will show that the 
algorithm is protected against the obvious sources of instabilities. This procedure can also be seen 
as an original application of the gradient flow, recently studied by Liischer [32| • 

This is not the whole story. The relative orientation between the canonical complex volume-form 
and the real volume-form, characterizing the tangent space of the thimble, contributes a phase to 
the integral. We can be sure that such phase is approximatively constant only where the quadratic 
approximation of the action is valid, but it is not obvious over the whole thimble. One expects 
such phase to change rather smoothly, and to affect only the sub-dominant part of the integral, 
but we have no clear evidence in support of this expectation. Moreover, the only procedure that 
we know to compute this phase scales very badly with the problem size. This difficulty definitely 
reduces the attractiveness of the algorithm that we describe in this paper. But, we believe that 
the approach that we have started to investigate deserves closer attention, since a better way to 
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deal with this residual phase does not look impossible. Moreover, the experience of the pioneering 
works [33| in lattice gauge theories suggests that very qualitative, but important, information on 
the phase structure might be gained already on tiny volumes. Unfortunately, the sign problem 
currently prevents the simulation even of lattices as small as 4 4 , in the high density regime. Any 
progress even on very small lattices might be very valuable. 

The paper is divided in two parts. In the first one (Sec. UTJ) we describe our approach in the case 
of a scalar field theory. In particular, we define the regularization in Sec. Ill Al and we analyze its 
symmetry properties and prove its perturbative equivalence to the standard formulation in Sec. Ill Bl 
In the same section we also introduce some basic Morse theory, in order to better understand the 
meaning of our formulation (Morse theory is not used as a justification, though). Finally, we 
illustrate and analyze the algorithm in Sec. Ill CI In the second part (Sec. IIIip we consider the 
extension to QCD^i. The definition of the functional integral needs to be adapted to the presence 
of local gauge invariance. The extension of the concept of Lefschetz thimbles to that case is standard 



3ll . 1 341 ] and it is presented in Sec. IIII Al The analysis of the symmetries and of the perturbative 



expansion is done in Sec. IIII BL where we also compare our approach to the standard one at = 0. 
Sec. IIIICI comments on the new aspect of the algorithm in the case of QCD/i. Finally, Sec. IIVI 
contains our conclusions. 



II. SCALAR FIELD THEORY WITH CHEMICAL POTENTIAL 



In this section, we consider the lattice discretization of a scalar QFT with chemical potential 
and quartic self interaction in d Euclidean dimensions. This theory represents a relativistic Bose 



gas at finite chemical potential and has a sign problem 35|, which can be successfully treated both 
via complex Langevin simulations 35], and with a reformulation of the path integral over current 
densities [36]. In this paper, we use it only as a simple framework to describe the fundamental 
features of our approach. 



A. Definition of the Path Integral 



In this section, we consider the model defined by the following lattice action: 

d-l 



S\ 



E 



(2d + m 2 ) <j>* x 4> x + A(C^) 2 - UZe-^fe+e + 



u=0 



(2) 



where <j> x is a complex scalar field, living on the sites x £ [0, L — l] d , and (ff x is its complex 
conjugate. The mass m, the coupling A and the chemical potential /i are all real and positive. A 
generic observable is computed as: 



{o)= \j c n^ e " 5W °[^ z =J c ii' 



(3) 



where the vector space C = C v ~ M. 2V is the domain of integration for the complex variables 
(f) x = < t >1 ' x ^ >2 ' x ; where the 4> a ^ x (a = 1,2) are real. Due to the presence of fi, the action S is 

not real, the real part of the Boltzmann weight 3i(e -S ) is not positive, and the system has a sign 
problem. Expressed in the real variables {<p a ,x}, the action ([2]) reads: 



S[Wa,x}] = £ 



' n n 
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d-l 



^ <t>a,x(f>a,x+v + ^2 i Smh |U E ah 4>a,x4>h^b ~ cosn A 4 ^ 



a i/=l 



a, 6 



(4) 



In order to introduce our approach, we first need to complexify the action ([2]). This is done 
by promoting to complex variables both the real part (4>i,x) and the imaginary part (4>2,x) of the 
field 4> x that enter in the formulation (jlj) (exactly as it is done in the case of the complex Langevin 
equation [351]): 

^^42+^2 (« =1,2). 

Inspection of Eq. shows that the action S[<p] is holomorphic in the (now complex) variables 
<fia,x, that parametrize the vector space C n , n = 2V. 

The equations of steepest descent (SD) for the real part of the action (Sr[((>] = $t(S[(J)])) are our 
second ingredient. They read: 



6<f>™ 
SSr[<P(t)} 



Va, x, 



Va, x, 



(5) 



Note that these are not the complex Langevin equations (at zero noise), which read: 



SS 



dr 



<3> 




+ 



SSr 

cj(S) ' 
S(paJ 

SS R 

aJ 1 ) ■ 



(6) 



Instead, the equations of SD can be reformulated, using complex variables, as: 

d 



dr 



5S[4>(t)} 

9a,x{ T ) = 7f ' Va ' X ' 



^<t>a,x{r) 



5S[<P(t)] 



Va, x, 



(7) 



For brevity, we also define the multi-index j = (R/I,a,x), which can be used to express the SD 
equations more concisely as <j>j = —djSft, and the multi-index k = (a, x), in which the SD equations 
become 4>k = —dkS. 

Finally, our approach (that will be justified only in Sec. Ill Bp consists in computing the observ- 
ables as: 



1 



0j Jo a 



S[4>]q\ 



Jo 



II- 



J a,x 



(8) 



where the set ^7o is an integration n-cycle defined as the union of all those curves that are solutions 
of the SD equations (0), or equivalently (|7j), and that end at the point = 0, in the limit r — > +oo. 

In the context of Lefschetz-Picard/Morse theory the point <ft = is called a critical 

point, because, by definition, it is a non-degenerate stationary point for the function S[4>], in the 
sense that 



6S[<f>] 



0, Va, x, 



and 



det 



J a,x |0 = o 



7^0. 



(9) 



y/ |0=o 
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The n-cycle J7b is the Lefschetz thimble mentioned in the title. One can prove [30] (Prop. 2.24) 
that Jq is a smooth manifold. Moreover, j7~o has real dimension n = 2V, as it should be in order to 
be an acceptable replacement of C. This can be seen as follows. The function S[4>] is holomorphic 
at <j) = 0, which is also a non-degenerate stationary point, because of condition ([9]). By Morse 
lemma (see e.g. 28], Prop 3.2, for holomorphic functions), there are holomorphic local coordinates 
<j){4>), where, in a neighborhood of <fi = 0, S[(f>] takes the form: 

s[4>)= E ^+- = E[(^ ) ) 2 -(tf ) ) 2 +^f ) tf ) )]+- 

k=(a,x) k 

The Hessian Hr of the real part of S[4>] can be seen as a matrix Hr G Hom(IR 2 " , R 2n ) (i.e. in the 
basis of the variables 0^ and 4> k ^), with precisely n positive and n negative eigenvalues. This 
implies that there are precisely n directions in which the SD flow comes from (f> = at r = — oo 
(the unstable thimble) and n directions in which the SD flow leads to (j) = at r = +oo, which 
defines our stable thimble J7b- Hence, the thimble J~q has real dimension n, precisely as the original 
integration cycle C. 

The fact that the imaginary part of the action, Si = 9(5*) = ^(S — S) is constant along the 
trajectories of SD is a straightforward consequence of Eq. (J7|): 

^ S R/i = \^ S± ^) = \ E (-9 k S-d k STd k S-d k S) = l -H^ll 2 

k=(a,x) 

Since S is continuous and well defined at = 0, it follows that Si is constant along the whole 
thimble J~q. This means that e tSl is an inessential constant phase, that can be factorized out of 
the functional integrals in Eq. (|8j). Moreover, Sr[4> = 0] is the absolute minimum of Sr along the 
whole thimble, and therefore the only non trivial part of the action is bounded from below. 

However, the fact that e lSl is constant, does not mean that the integrand is real and positive. 
In fact, the measure Y\ ax d(pa,x stands for the complex canonical volume- form in C n , that needs 
to be evaluated on a basis of the tangent space T^Jq) of J$ in <fi. Since T<f,(Jo) is a real n- 
dimensional vector space C C n , its basis is not necessarily aligned with the canonical basis of C n . 
This misalignment produces a phase, which is univocally determined by T^Jq). In Sec. Ill CI we 
will describe a (rather expensive) procedure to calculate it numerically. 

Now that we have defined precisely what we want to compute, we need to address two obvious 
issues: the first one is how to justify the study of Eq. (|8|). In fact, the integral Z in Eq. ([3]) does 
not coincide with Zq, defined in Eq. ([8]). Nevertheless, we will argue in Sec. Ill Bl that the system 
in ([1]) is physically at least as interesting as the original formulation. The second issue is to find 
an algorithm to compute Eq. (JH]). In fact, it is far from obvious how to perform an importance 
sampling of the field configurations in J7b- An algorithm is proposed and analyzed in Section III CI 



B. Justification of the Approach 

The main justification to study the system defined by Eq. ([8]) is the fact that the latter defines a 
loca0 QFT with exactly the same symmetries, the same number of degrees of freedom — belonging 
to the same representations of the symmetry groups — and the same local interactions as the original 
theory. Moreover, at jj, = the cycle Jq coincides with C. Furthermore, the perturbative expansions 

1 One might worry that the definition of the integration cycle introduces a subtle non-local interaction, since the 
allowed values for the field <f>k at one space-time point k depend on the values of the field at the other points. 
But the thimble Jo is just one representative of a homology class. As long as the complementary of the set of 
critical points is open, it is certainly possible to deform Jo (without changing its homology class) into a new cycle 
that — in the neighborhood of any given configuration (f> g Jo — has the form ®% =1 Vk, where <f>k £ 14 C Ck Vfc, thus 
removing the suspicion of non-locality. 
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of the two systems coincide at every /i. Since universality is not a theorem, we certainly cannot 
prove that Zq and Z are physically equivalent. But given the above properties, it would be 
extremely interesting if they did not, as it would show that those properties are not sufficient to 
characterize a QFT uniquely. 

In the following Sec. Ill B II we analyze the symmetries of the new formulation. In Sec. Ill B~2l we 
show the perturbative equivalence of the two approaches. Then, in Sec. Ill B 31 we introduce some 
rudimentary Morse theory. This will not enable us to determine the exact relation between the 
integral Z over C and Zq over j7o, but at least it represents a convenient framework to gain insight 
into the relations that one might expect between the two formulations. 



1. Global Phase Symmetry 

The only symmetry of the action in Eq. ([2]) that could be non-trivially affected by the sub- 
stitution of the integration cycle C with j7o, is the U(l) symmetry associated to a global phase 
rotation: (ft x — > e ta (j) x . This symmetry holds for any \x. In the complexified system, such symmetry 
translates into an 50(2, C) symmetry associated to the rotations of the new (complex) variables 

\x \ . .n<r, / "'!..,- 



e aa2 Y' x . (10) 
n,x j \ <t>2,x ) 

These transformations leave the point eft = invariant, and hence do not produce zero eigenvalues 
in the corresponding Hessian matrix, which is consistent with the observations made above. 

When considering the path integral defined over j7o, we need to ask whether the physical U(l) 
symmetry is preserved. In the complexified action ([2j), the original U (1) symmetry group is mapped 
into the real subgroup 50(2,M) of the whole 50(2,C), which corresponds to real values of 0.02- 
The action ([2]), which is invariant under 50(2,C), is obviously invariant also under its subgroup 
50(2, R), but it is not obvious whether any of these symmetries is defined in j7o, i.e., whether the 
configuration <ft = e ar72 (ft belongs to Jq, whenever (ft does. 

By definition, (ft 6 Jq if there exists a curve (ft{r) that solves the SD equations (|5l7p . with 
(ft(0^ = (ft and <ft(+oo) = 0. If so, we can define the curve (ft(r) = e a<J2 0(r), which obviously starts 
at (f)(0) = (ft and ends at (ft(+co) = 0, and it also solves the SD equations ([7|), in fact 



dr dr dej) 

In the third step we have used the reality of aa2, for rotations in 50(2,R), while in the last step 
we have used the covariance of the gradient of 5 under the transformation (jlOp . By uniqueness of 
the solution, we conclude that <j) £ Jo whenever (ft does. Hence, the formulation based on J7b has 
the same symmetries of the original formulation. 



2. Perturbative Analysis 



We consider now the perturbative expansion of the system defined by Eq. ([8]). In a nutshell, the 
reason why the perturbative expansions of Eq. ([8|) and Eq. ([3]) coincide is the fact that Gaussian 
integrals (times polynomial P) in the complex plane 



/ dz e~ z2 P(z) 

7 
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are independent on the path 7, as long as 7 joins the region at infinity where | argz| < ir/4 with 
the other region at infinity where | argz| > 37r/4H 

In more detail, when expanding in powers of A the expression ([8]), we have to take into account 
the fact that Jq = Jo(\, //) depends itself on A. Hence, at perturbative order p, we have to consider 
expressions like: 




S[<t>;\,iA 




|A=0 



for a generic observable O. This expression generates terms like 

, , dP 



Mo,n) d\P\\=o 



-S[4>;\,p]q 



and genuinely new terms like 

d 

d\ 



(11) 



(12) 



(13) 



where P is some polynomial in eft, whose coefficients depend on [i. The terms like Eq. f)12[> are 
Gaussian integrals (times polynomial) defined over an integration cycle J7o(0,//), which is defined 
as the path of SD for the free part of the action at finite fi. Since for Gaussian integrals the non 
trivial element of the homology class is unique, the integral (I12D coincides with the integral of 
the same function along C, assuming that the latter converges, which is true as long as standard 
perturbation theory is well defined. 

Consider now the terms like the one in Eq. (|13p . There, the only dependence on A appears in 
the integration cycle. Hence, Eq. (113j) measures the variation of the integral under infinitesimal 
variation of the integration cycle around j7o(0,/i). But the cycle J7o(0>/-0 corresponds to the path 
of steepest descent for the integrand e~ s ^'°>^ times polynomials. In particular, the cycle J7b(0, fi) 
lies in the interior of the region of convergence, and the integral that is differentiated in (|13|) cannot 
change for infinitesimal variations, which means that contributions like f)13[) always vanish. We 
have proved that the perturbative expansions of Eq. (|8j) and Eq. ([3]) coincide. 



3. Some Insight from Morse Theory 

It is usually very difficult to make any analytic statement about a QFT beyond symmetries, 
locality and perturbation theory. These are actually the properties that usually justify a legitimate 
regularization of a QFT. Nevertheless, it is interesting to investigate further the relation between 



the path integral defined over C and the one over Jq. In fact, Morse theory [30|, [3l( can be used 
to gain much insight into this question, although we won't be able to establish any exact relation. 
In this section, we summarize some general results from Morse theory (see, in particular, Sec. 3.2 
of [lH), with a special attention to the cases of interest for us. 

Let S{x) be a complex function of n real variables (27, . . . , x n ), that we analytically continue to 
complex values — > Z/.. Assume that S(z) has only finitely many critical points, and it is generic 
enough that they are all non-degenerat^l (i.e. it is a Morse function). As already mentioned, these 
are points where the gradient of S vanishes, but the determinant of the Hessian is non zero. For 



In more homological terms, such paths define the only non-trivial element of the relative homology class 
Hi(C, Ct;Z), for a sufficiently large T. The set Ct is the set of all z G C such that 5R(z 2 ) > T. For a defi- 
nition of relative homology, see, e.g. [37t ] . 
3 Later we will consider sets of degenerate critical points, which typically appear in presence of symmetries. 
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each critical point z = (f) a , a G E, the Hessian H a of Sr = ?R.(S) at (j) a can be seen as a bilinear 
real form in the In- variables (uk,Vk), where Zk = Uk + ivk- Hence H a has precisely n positive and 
n negative eigenvalues. 

To each critical point <j) a we attach a stable Lefschetz thimble J a , defined as the union of 
all flows that satisfy Eq. ([5]) and tend to </> CT when r — > +oo. For each a, we also introduce 
an unstable Lefschetz thimble /C CT , defined as the union of all flows that satisfy Eq. ([5]) and go 
to 4> a when r — > — oo. Because of the holomorphicity of S in <p a , all such thimbles have real 
dimension n. Moreover, for a generic choice of the parameters in S, all J a and JC a extend to 
infinity without crossing other critical points. When this is the case, the J a provide a basis of the 
relative homology group := H n (C n ,(C n )T;%), while the K, a generate the relative homology 
group H~ := H n (C n , (C n )~ T ;Z)@. There is a duality [27J between the group and the group 
H~ , which is realized by the bilinear form: 

( , ):H+®H~^Z, (14) 

which associates to each pair of cycles C G and C" G H~ the (oriented) intersection number 
of the two cycles. In fact, in the generic case, two half-dimensional manifolds intersect in a zero 
dimensional set. In particular, the basis J a and K, a are dual to each other, because their intersection 
is either zero or amounts to the single point </> CT , and hence {J p ,fC a ) = 5 Pj(T . These observations 
lead to a general formula that enables us to express a generic integration cycle C G in the basis 

{JaY 

C = Y,KaJa, (15) 

where n a = (C,/C CT ). This formula implies that, in order to reproduce exactly the original integral 
over C, we should consider not only the cycle j7o, but also the contribution (with sign) from all 
the critical points of S[<fi] in C". However, the argument of [311 ] suggests that most of these other 
critical points might give either an exponentially suppressed contribution or no contribution at all. 

The argument goes as follows. Let s mm = min^ g c Sr(4>; /i) be the global minimum value of the 
real part of the action in the original manifold. The full set of critical points £ is the union of the 
three disjoint subsets: 

S = {a G S | <P„ G C} 

£< = {a G S | CT i C & SflOM < s min } 

£> = {a G S | CT i C & S^) > s min } 

The critical points in <p a G S< do not contribute to Eq. (fT5|) . because in (p a the value of 
Sr is already lower or equal to its absolute minimum in C, and Sr can only decrease further in 
the unstable thimble K a . Therefore, ¥i a can never intersect C and n a = 0. The critical points 
^eS> may or may not contribute, but their contribution is exponentially suppressed by a factor 
e ~s R (4>^)+s min ^ w fth respect to the thimble associated to the absolute minimum of Sr in C. One 
expects, generically, this suppression to be further enhanced toward the infinite volume limit, but 
we cannot exclude that, for example, a large number of critical points might form and accumulate 
near the absolute minimum in C. 

Finally, the critical points 4> a G So necessarily contribute with n a = 1, because the thimble /Co- 
intersects C precisely once in <p a , but their contribution is suppressed, if they are not global minima 
(for the same argument used for points in £>). Actually, in the case of the action ([2]), there are no 



4 The space Xt (resp. X T ) is defined as the set of those points of X such that > T (resp. Sr < —T). 
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further stationary points in C, besides = 0, but there must be many in QCD (sec Sec. IIII B 3|) , 
This argument is not conclusive, since we cannot exclude, for example, an accumulation of critical 
points near the global minimum in C. However, it shows that assuming a regularization defined 
only on the thimble associated with the global minimum in C is not in contradiction with anything 
we know from Morse theory. 

It is also interesting to consider more closely the case of a QFT with chemical potential ji. 
In particular, it is interesting to see what happens when we start from a real theory and switch 
on jj,. Consider first the case of fj, = 0. In this case, the action S[<p] is real, and the condition 
of stationarity SS 1 ^] = imposes n equations with n unknowns. Hence, we expect, in general, a 
discrete set of solutions. These stationary points can be minima (local or global), maxima or saddle 
points. Notice that at fj, = 0, the flow of SD preserves C. If the system has just one minimum (f)Q, 
and no saddle points, the manifold C coincides with the thimble J fa. The saddle points that might 
be in C represent stable limits only for a set of zero measure in C; all the other points will eventually 
flow to (fro. Hence, even in the presence of saddle points, the closure of J fa still coincides with 
C. On the other hand, if S[(j)] has further (local) minima in C (besides ^o) the situation changes. 
Each of these minima represent the stable limit for a measurable subset of C. The contribution of 
these subsets to the integral is exponentially suppressed, for the same reason explained above, but 
it might be important, if they are many. In summary, at [i = 0, we can write: 

C = 3v modulo a set of zero measure. 

<J | (f>(j is local minimum for S 

This shows that, even those critical points $ a where n a ^ may actually give a vanishing contri- 
bution to the integral. 

When we switch on /i ^ 0, the action S[(j)] becomes complex, and the condition c^S 1 ^] = 
in C becomes a system of 2n equations with n unknowns. As a result, all the generic stationary 
points in C are shifted outside C, unless some symmetry protects them. Does this mean that all 
the minima that contributed at fj, = suddenly become irrelevant as soon as \i is switched on? If 
not, do all the stationary points in the nearby of C suddenly become equally important? Morse 
theory enables us to verify that the transition is actually smooth. 

In order to illustrate the mechanism of this transition, assume that the following ansatz is valid 
for small values of //: 

where S°[(p] and S [4>] are real functions of the variables (j) = {4>k}k=(a,x)-> when the cfi are real, and 
holomorphic for complex values of <fi. 

An expansion to first order in [i shows that, if (j) is a critical point at [i = 0, the new critical 
point (ft = (ft + 5(j) is shifted by: 



ill [diS 1 ^]) (d 2 S°[4>. 



-l 

Lk 



In the variables cj) = <fr( R ' + i<p^> this reads: 



J k 



-1 

' l,k ' 

In order to compute the contribution of J? to Eq. (fT5j) . we need to compute the index n^. For 
this we need to check whether the unstable thimble /Ct, associated to (ft, intersects the original real 
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manifold C, which is defined by (ft^ = 0. In the neighborhood of (ft, the unstable thimble 



t— oo 



which is valid only for t, \ \x\\, \\y\\ up to 0(/j,). Now we can appreciate the different fate of the local 
minima in C from the saddle points. If (ft is a minimum at fj, = 0, the eigenvalues of d 2 S°[(ft] are 
positive and IC-g is characterized by X = 0. Hence, there is always a choice of Y that cancels the 

shift 5(ft^ and /C^ intersects C. If instead is a saddle point at /z = 0, the possible Y are restricted 

to the eigenvectors of d 2 S°[(ft} associated to positive eigenvalues. Hence /Ct does not intersect C 
for a generic choice of the parameters. This is compatible with the expectation that, generically, 
no dramatic change happens at small fi. In fact, what gives a finite contribution at fi = (the 
local minima) still contributes at small /i, since ^ 0. While the terms that have zero measure 
at n = (the saddle points), have, generically, also = at small \i. 

This last discussion is obviously relevant only to show the consistency of the picture at small 
fx, and certainly not to justify our approach at finite fi. The justification of the latter relies on the 
considerations done previously in this section. 



C. Algorithm 

In this section we describe an importance sampling Monte Carlo algorithm to simulate the 
integral in Eq. ([8]). In this section we ignore the phase due to the measure, which will be taken 
into account in Sec. IIICTI via a reweighting step. 

Since the imaginary part of the action Sj is constant along J7o, we can rewrite Eq. ([8]) as: 

[ T\dcft x e~ s ^O[(ft} = ^e-^ f T\d(ft x e~ s ^O[(ftl (16) 
z J Jo V A) J J-, V 

and the phase factor e~ lSl effectively cancels from the expectation values. Hence, we need an 
algorithm to simulate the real action Sr on Jq. Note that Sr is bounded from below on J§. 

We would like to compute the integral (|16p through a Langevin algorithm, constrained in J§. 
The corresponding Langevin equations are: 

d AR) _ _ SS R , (R) 

dT 9a,x — —m ' r la,x, 

D( Pa,x 

which we summarize hereafter as <ftj = —djSn + m, where j = (R/I,a,x) is the multi-index 
introduced earlier, and rjj is a random field with the usual properties: (r/j) = 0, (rjjrjji) = 28jj>- As 
already noted, these are not the complex Langevin equations. Instead, (fT?]) coincides (for rj = 0) 
with the equations of steepest descent for Sr. As opposed to Eq. ©, which might develop non- 
trivial attractors in C n , Eq. (|17[) would drive the system all the way down to Sr — > — oo, unless it 
is constrained in J$. 

We now come to the problem of constraining the system in Jq. The drift term in Eq. (|17h 
keeps the configuration in Jq by definition. The difficulty lies in extracting a noise n tangent to 
j7o, despite the fact that we lack a practical local characterization of J§. In fact, if we start from 
a configuration (ft € Soi it is not clear how to determine which directions from (ft are tangent to Jq 
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and which are orthogonal to it. This depends on the long time evolution of the nearby paths of 
steepest descent, which is difficult to determine locally. 

On the other hand, the space T^ =0 (j7o), tangent to Jo at (ft = 0, is easy to compute. Hence, 
the random noise can be projected onto the correct subspace Tq(J'q), at (ft = 0. After that, it 
can be parallel transported along the flow dS R , that connects = to a previously generated 
configuration (ft, and it can be added to (ft. The combination of the Langevin noise steps with drift 
steps drives the system naturally toward the regions of Jq that dominate the functional integral. 

The concept of Lie derivative provides the natural tool to parallel transport a vector r\ along 
the flow dSft. As we shall see below, this is also straightforward to implement. In this way, the 
importance sampling is realized in the usual sense of the Langevin equation, that relies on the 
correct balance between a drift term and a noise term. There are important questions of stability 
in this procedure, and we address them in detail in Sec. Ill C 21 

The following algorithm is, essentially, a Langevin algorithm except that each time we want to 
add a noise vector, we have to transport it back and forth to the origin, in order to ensure that it 
belongs to T^Jq). The detailed procedure is the following (steps 1-7 are preparatory; steps 8-14 
should be iterated): 

1. Compute the Hessian matrix (d 2 S R ) := (d 2 S R )^ =0 at the critical point (this can be done 
analytically once and for all, and it is reported in Appendix I A 1|) . 

2. Extract a random field rj (with 2n real components), from an isotropic distribution (normal- 
ization will be done later). 

3. Project the noise rj into the eigenspace of (d 2 S R ) of real dimension n associated to positive 
eigenvalues of (d 2 S R ) @. In this way, we obtain a vector 7711 , parallel to Jq in the origin (ft = 0. 

4. Normalize the vector 7711 such that ||r?|||| = e. The sphere of radius e must be sufficiently 
small so that the second derivative of the action can be approximated by a quadratic form 
in the fields (ft. This is the region where the Hessian matrix computed in (ft = can be used 
reliably. 

5. Evolve the vector 7711 with the equations of steepest ascent. For this first step, it is irrelevant 
how long we follow the curve (say for time To): this produces the starting configuration <fto. 

6. Extract a new random field as in step [2j 

7. Transport r/ 1 ) along the path of SD that brings the configuration (fto back to the sphere of 
radius e. This is done by ensuring that Lie derivative CQs R (r/^ (r)) of r/ 1 ) along the flow 
defined by dS R is zero. In fact, the condition that the Lie derivative is zero is equivalent 
to ensure that the transported vectors rj^^r) commute with the field dS R , and hence, the 
paths that follow the two vector flows in different order commute. Evolving (t) while 
ensuring Cgs R (ri^ (r)) = is also straightforward to implement numerically: 



= W? (1) (^)) = 

= [dS R , V W(T)] = 

= Y J d 3 S R d 3 r ] f{r) -Y^rif\r)d j d j ,S Rl 

j j 



5 We remind the reader that (d 2 S R ) has n positive and n negative eigenvalues, by holomorphicity. 
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which is equivalent to: 

^vf(r) = ^24\r)d r d j S R , (18) 

j' 

and can be solved numerically by applying standard methods of integration of ordinary 
differential equations (ODEs). Note that (t) remains constant (in direction) only if it 
happens to be an eigenvector of 8 2 Sr. 

8. Project rf 1 ^ onto the positive eigenspace of (d 2 S R ) to produce rfjp. 

9. Transport r/j; along the SD curve that leads from the sphere of radius e to 0o- This is done 
by ensuring that the Lie derivative of 7]^ along the field 8Sr remains zero, as described in 
step [71 Note that this means, in particular, that (r) remains tangent to J§. 

10. Once rh (t) has been evolved up to To, it can be added to <pQ. The norm of r/y (r) is 
determined by the theorem of stochastic quantization: it must be sampled from a suitable 
distribution (e.g. Gaussian) with standard deviation equal to \2d~t. 

11. Perform one Langevin (i.e. steepest descent) step of length dt. This produces a new config- 
uration (jP-i. 

12. Extract a new noise r/ 2 ) . 

13. Evolve via SD down to the origin and check that it meets the ball of radius e and that 
it falls into the positive eigenspace of (d 2 Sj i ). If not, reduce dt and repeat. 

14. At the same time transport t/ 2 ** along the path connecting <f>^ to the origin, as described 
in step [7J @ 

15. Iterate from [8l to 1141 ad libitum. 

1. The Residual Phase 

The algorithm described above only samples the configurations; it does not yet take into account 
the phase that comes from the misalignment of the tangent space T^Jq) with respect to the 
canonical complex basis, in which the complexified integral is formulated. In order to do so, we 
need to compute an orthonormal basis of TA^Jq), for each configuration eft that we sample, in terms 
of the canonical basis and compute its determinant. As already noted, the tangent space is easy to 
compute only in = 0, and it can be computed in other configurations only through the parallel 
transport along the flow. Unfortunately, in this case, Liouville's formula cannot be used directly to 
transport the determinant along the flow, and we currently see no better option than transporting 
every single vector of a full basis, using Eq. (fT8|) . Such procedure costs 0{V 2 L^) both in terms of 
storage and flops, where V is the four dimensional volume, and L5 is the number of steps in which 
the SD flow is discretized. Moreover, the calculation of the determinant requires 0(V S ) flops. 

6 One could imagine generating random noise vectors directly at the origin <f> — and transporting them only 
upwards along the directions of steepest ascent, thus omitting steps [7| and 1141 However, isotropy of the noise 
would not be guaranteed. This might not be a problem for some algorithms that do not require isotropy of the 
proposal for correctness, but only detailed balance. This might be the case if we add an accept /reject step as in 
the Langevin Monte Carlo algorithm ; 38, 39]. However, it is not completely clear to us, whether detailed balance 
is satisfied in this case. So, we require isotropy of the noise in this paper, which is ensured by the above procedure 
and is compatible both with a Langevin and a Langevin Monte Carlo algorithm. 
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This is of course a lot. One could argue that this is already much better than the 0(e ), which 
is expected from a direct simulation of the model. But this clearly limits the approach to very 
small lattices. The importance of simulating at least small lattices should not be underestimated, 
however. In fact, the real trouble with the sign problem is not only the bad scaling, but also the fact 
that even lattices as tiny as 4 4 appear intractable by brute force. This is unfortunate, because the 
experience from the early days of lattice field theories, suggests that very qualitative, but useful, 
information on the phase structure might be gained already from tiny lattices [33]. It becomes 
then crucial to understand how much the sign problem is reduced with the method proposed here. 
To this purpose, we can presently only argue that the phase must be essentially constant over the 
portion of phase space that dominates the integral, and the fluctuations should only determine 
the corrections to the dominant behavior. However, we are unable to provide clear quantitative 
support to this qualitative argument. This question should be definitely assessed through tests. 
On the other hand, we observe that this approach is quite new in many respects, and we expect 
that new ideas might solve or substantially improve on the problem described in this section, even 
before proceeding to expensive tests. 



2. Remarks on Numerical Stability 

The algorithm described in the previous section requires the solution of a few systems of ODEs 
and it is mandatory to comment on their expected numerical stability. 

Step [5] consists in integrating the equations of steepest ascent evolution from the sphere of radius 
e towards the interior of J§. Such integration should be stable against perturbations, because the 
evolution in those directions should suppress the eigenmodes leading out of Jq (this is exactly true 
in the quadratic approximation of the Hessian). Moreover, step [5] needs to be performed only once 
in the initialization phase. Steps El M and HH require the solution of ODEs along a path in j7o 
which is already known. The components of the vector n which are orthogonal to J7b are likely 
to be enhanced when approaching <p = in the descending direction, but these components are 
projected out when the point <ft = is reached within the precision e. Certainly, the ODE in 
Eq. (fl~8j) is expected to be stiff, and the integration method must be chosen accordingly. 

The integration in step [13] deserves more attention. In fact, when we try to construct a new 
path that follows the curve of steepest descent that goes from a point in the interior of j7o towards 
the point = 0, any small perturbation outside J$ is expected to be associated to diverging 
eigenmodes that drive the system away from = 0. There is, however, considerable experience 
in dealing properly with these kinds of ODEs which are boundary value problems rather than the 
more common initial value problems. In these cases it is essential [13] to use the information at 
r = oo. Once this constraint is imposed, the ODE can be solved via a finite difference method, 
which involves exactly the same kinds of derivatives that are typical of the lattice discretization of 
QFTs. In order to have a chance of success, it is crucial to start from a good initial guess. This 
is available in our case. In fact, from the previous path and from the noise vector projected along 
T(f>{t)(3o) i w e can propose a first guess not only of the new configuration 0(t + dt), but also a guess 
for the whole new path that joins 4>(t + dt) to = 0. This is done by exploiting the fact that the 
noise field defined along the SD path (as constructed in step [9|) commutes with the SD flow, and 
hence we can add to each point of the previous path the corresponding vector rj, which produces 
a guess for the new path which is good to 0(dt 2 ). A final check is to verify that the imaginary 
part of the action is indeed constant. In other words, the integration path obtained in the previous 
step provides a very good guess for the subsequent path, and the solution of the boundary value 
problems prevents the accumulation of errors. 

In conclusion, the algorithm outlined above includes two tunable parameters: dt and e. The 
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step-size dt needs to be sufficiently small to ensure correctness of the algorithm and sufficiently 
large to enable an effective sampling of the fields. On the other hand, the ball radius e needs to be 
sufficiently small, so that the quadratic approximation of the Hessian is valid, but also sufficiently 
large, such that the small 0(dt 2 ) perturbations to the paths in J7o do not cause the paths to miss 
the sphere around the origin. Only suitable numerical tests on realistic conditions can tell whether 
such compromises are possible. 



3. Costs Estimate and Possible Improvements 



The cost of computing the determinant of the tangent space T^Jq) has been already estimated 
in Sec. Ill C II This is by far the dominant cost. However, in the optimistic perspective that the 
computation of such phase may be simplified or avoided, it is interesting to estimate also the other 
(presently sub-dominant) costs. 

The next dominant part of the algorithm is the solution of an ODE of a system of size V for a 
length L5. If it is carried out with a finite difference method, this task costs 0{VL^) in memory 
and 0(V X571FD) in flops, wher e ny r> is the number of iterations used by the solver. Further insight 
is provided by the observation [32] that the distances in the fifth direction have physical dimension 
[length] 2 . Hence we might expect the unfavorable scaling L5 ~ y 2 l d . On the other hand, the 
fifth dimension does not entail a true quantum dynamics, as the fields at finite r are completely 
determined by those at r = 0. For this reason, the autocorrelation length cannot be affected by 
this growth of the problem. Moreover, it is very likely that methods of over-relaxation 4l|, |42|] and 
Fourier acceleration 43J will be important to integrate the equations of steepest descent /ascent 
efficiently. Further possibilities to accelerate the evolution in the fifth dimension will be mentioned 
in the case of QCD. 

In any case, the algorithm proposed in this paper contains a number of new elements with 
respect to the well known and reliable tools to which lattice QCD theorists are well accustomed. 
In this context, one can imagine many unexpected difficulties, but one should also expect that 
other computational sciences will suggest new strategies to overcome them. In any case, it is easy 
to foresee a hard work of testing and tuning. 



III. QCD WITH BARYONIC CHEMICAL POTENTIAL 

In this section, we apply the same analysis done in Sec. [II] to the case of QCD^. 



A. Definition of the Path Integral 

1. The Standard Formulation 



For the lattice action of QCD at finite density, we assume the classic Wilson regularization 
4J, |45J] and introduce the baryonic chemical potential as usual 461 ] : 



S[U] =pJ2 



1 _ Isft Tr U^ x ) 



+ N f Ti log Q[U], 



(19) 
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where, Nj is the number of degenerate quark flavors, U^ jU (z) is the classic Wilson plaquette, and 
the Dirac operator Q is defined as: 

Q[U) xy = (m + 4r)d xy -± £ (r-^+^e^U^x)- X - £ (r + lv )5 x ^ y e-^XJ v (x-P)-\ 

^=0,3 ^=0,3 

(20) 

where we assume either periodic or anti-periodic boundary conditions and a finite box of volume 
V = 7 3 x T. In the standard formulation, the observables are defined through the path integral: 




where C = {U u {x) £ 517(3), Vi 6 [0 . . . L - l] 3 x [0 . . . T - 1], v = . . . 3} is a manifold of real 
dimension n = V x d x (Af? — 1) = V x 32. Note that we have explicitly performed the integration 
over the fermionic quark fields, to produce the exact effective action in Eq. (|19p . Finally, note that 
the gauge does not need to be fixed, although it is possible to do it. 



2. The New Formulation 



In order to introduce our formulation we need to define the corresponding integration cycle J7b- 
The first step is to complexify the system. This is achieved by extending the SU (3) gauge group to 
SL(3, C). This correspond^!] to the complexification of the algebra su(3) 4xV — > st(3, C) 4x v , which 
is a vector space of complex dimension n = 32 x V: 

A a u (x) -> A a u > R (x) + iA a /{x) a = 1 



N 



1. 



In Sec. IIII B 41 we discuss the domain of holomorphicity of S[A] in greater detail. Here, we simply 
observe that S[A] is holomorphic, as a function of A, in a neighborhood of A = (except for a 
discrete set of values of n that will be discussed later). 

The second ingredient that we need is a notion of derivative with respect to the fields U v {x) 6 
SL(3, C). The natural definition of field derivative on a lattice regularization of QCD is the left 
covariant one. In order to derive the analogues of Eqs. ([5]) and ((ZJ), it is convenient to define the 
following set of left covariant derivatives: 



d 



V^ a F\U) :=-F[e^U u (x)] la=0 . 
V x ,u, a F[tf] := 0, 







V*„, a F[E7] :=— F[^U v {x)] 



V x ,u, a F[U] := 0, 



da 

V R F\U ] ]:=—F 
Vl,v, a F[U] :=^-F [e- aT «U u (x)] 



|a=0 ' 

U u (x)^e- iaT ' 



(22) 



WP7t] := If 



U v (xye 



\ e -iocT a 



|a=0 



da 

V^tU] = ^F 



u v (xy e 



t e -aT a 



|a=0 
|a=0 ' 

|a=0 ' 



where T a are the (Hermitian) generators of the algebra su(3) in the fundamental representation 
(normalized as Tr(T a T{,) = |<5 a ,&)- Note that the derivatives above satisfy the Cauchy-Riemann 
relations on the functions that depends only on U or W , and that|f|: 



(23) 



7 Also in this case, the procedure of complexification coincides with the one adopted in the context of complex 
Langevin formalism [121 ]. 

8 We use the multi-index k — (x, u, a). 



17 



The important advantage of these definitions is that the derivatives of gauge invariant functionals 
are exactly gauge covariant — even at finite lattice spacing a. One should keep in mind that the 
derivatives (|22p do not commute. Instead, they obey the following commutation relations (which 
hold for any of the above derivatives): 



y,<r,t 



0~x,y3u,a fabc^ x,v,a 



where the f a b c are the structure constants, defined by: [T a ,jy = ifabcXc- Note, however, that the 
Hessian matrix of a function F[U] is well defined and symmetric at any of its stationary points. 

We will also need to express the derivative of a function F : <SX(3, C) n — > C along a 1-dimensional 
curve U(t) C SL(3, C) n (we suppress inessential (x, v) indices in this discussion). To that purpose, 
note that if U(t) is such a curve, generated at r by an infinitesimal left-translation in the direction 
of a a T a , i.e.: 



U(t + dr) = e 



dr a a iT a 



U(r) 



we can express: 



Hence, 



O'n 



-2Tr 



dr 



-F[C/(r)] = V a F[U(r)} 



-2Tr 



iT a {-U{r))U(r)- 1 



(24) 



We now need to discuss how the appearance of local gauge invariance, in Eq. (|19p . affects the 
construction of the manifold Jq. In fact, the point A = (as well as any other stationary point) 
changes non-trivially under general gauge transformations, and hence it cannot be an isolated 
stationary point. More precisely, every stationary point of S[A] belongs to a manifold of stationary 
points and, in particular, the Hessian is degenerate. In this case, the concept of (un)stable thimbles 
attached to A becomes ambiguous. 



The appropriate way to deal with these cases is explained in [311 . |34| : in presence of symmetries 
that act non-trivially on critical points, it is convenient to generalize the concept of a non-degenerate 
critical point into that of a non- degenerate critical manifold [13]. A manifold N C C is a non- 
degenerate critical sub-manifold of C for the function F : C — > M if: 

1. dF = along Af; 

2. The Hessian d 2 F is non-degenerate on the normal bundle v{N). 

Under these conditions, we can decompose the bundle normal to M as v{M) = v~{N) ® v + {M), 
where the first (second) bundle in the sum is associated to strictly negative (positive) eigenvalues 
of d 2 F. 

In the case of the QCD lattice action in Eq. (fT9|) . the manifold M represents a full gauge orbit of 
stationary points^!, which has real dimension uq := dim^Af = (V — 1) x (AT? — 1). As we complexify 
the system, the manifold N also extends to a larger manifold Afc of complex dimension uq. The 
manifold Mc is the orbit generated by application of all possible SX(3, C) gauge transformations to 



In the pure Yang- Mills case, JV also includes toronic degrees of freedom [is|, but this degeneracy is removed by 
the fermionic part of the action, as confirmed by the computation of the Hessian in the Appendix. 
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the configuration A = 0. Hence, the Hessian of the real part of a holomorphic and gauge invariant 
function F : C — > C can be regarded as a real matrix in Hom(R 2n , M 2n ), which has n — uq positive, 
n — no negative and 2uq zero eigenvalues. As stressed in [31] (see, in particular its Sec. 3.3), 
the n-dimensional integration cycle that we need should be build out of the stable manifold of 
curves of SD attached to a middle- dimensional manifold contained in J\fc- A natural choice for the 
middle-dimensional manifold in Mc is the original N itself. 

Finally, we need are suitable SD equations. The generalization of Eq. ([7]) to the left-covariant 
case leads to: 

^U u (x; t) = (-iT a V x ^ a S\U])Uu(x; r) (25) 



Similarly to the scalar model, Eqs. (J25J) are equivalent to minimizing the real part of the action 
Sr[U]. Moreover, the imaginary part Si[U] is conserved along those curves. Both of these proper- 
ties can be verified by using Eqs. ([HP and ([25]) : 

d „ 1 d ,„ , ^ 1 \ - , „ „ = „ _ „ . f - II VS II 2 



J? Sr/i = \^ {s±s) = \ S {-^kS-V k S^V k S-V k S) = i 



k=(x,i/,a) 



o 



After this long preamble, we can finally define the integration cycle J7b as: 

Jo : = \u € (5L(3,C)) 4y | 3C7(r) solution of Eq. (Ell) | 17(0) = U & lim f/(r) G AA (0) ) , 

(26) 

where J\f^ is the critical manifold that contains the point A = 0. The definition (|26p ensures that 
is an integration cycle of the right dimension (n — no) + (2wg)/2 = n. Moreover, the choice of 

the critical manifold AA°) ensures (as shown in Sec. IIII B 2|) that the perturbative expansion of the 

new formulation coincides with the standard one. 

Substituting C with Jq in Eq. (f2T|) concludes the definition of our procedural: 

(O)o = ^r[ \[dU v (x) e- s MO[U], Z = f ]JdU v (x) e~ S W, (27) 

In the next sections we justify why this new formulation is physically relevant, and propose a Monte 
Carlo algorithm to study it numerically. 



B. Justification of the Approach 

As already explained in Sec. Ill Bl we do not attempt to derive an exact relation between the path 
integral on the cycle C and the one on the cycle J7o- Our motivation to study QCD on the thimble 
j7o is that it is a non-perturbative definition of a local QFT with the same algebra of operators, 
the same degrees of freedom, the same symmetries, and the same perturbative expansion as QCD. 
If the continuum spectrum of QCD is an unambiguous prediction of these properties — as it is 
generally expected on the basis of universality — then studying the formulation in J§ is physically 
very significant. If that should not be the case, it would represent a very interesting surprise, and 
a major step forward in our understanding of QFTs. 

Motivated by these ideas, we examine, in the following sections, the symmetry properties 
(Sec. IIII B 1|) and the perturbative expansion (Sec. IIII B 2|) of our formulation. In Sec. IIII B 31 
we define a strategy to compare precisely the formulations in C and Jq at fi = 0. Finally, in 
Sec. IIII B 41 we comment on the branches of the logarithm that appear in the fermionic effective 
action. 

10 The measure dU v {x) needs to be evaluated on a basis of the tangent space of j7o, which produces a phase, as 
discussed in the scalar case. 
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1. Gauge Symmetry 

The only new symmetry that deserves special comments, in the case of QCD, is the SU{3) 
gauge symmetry. The SD Eq. (|25p is exactly covariant under gauge transformations U v {x) — > 
K{x)U v {x)K{x + v)~ l , but only if the transformations A belong to the SU{3) sub-group of the whole 
SL(3, C) symmetry group that emerged after complexification. This is due to the fact that — as 
opposed to the complex Langevin equation — in the SD Eq. (|25j) the conjugate term (T a V x ,u,aS[U]) 
appears, which transforms as: 

S[U}) (A(x)~ i y (T a V x ,u, a S[U))A(x)l 

More precisely, we can show the SU(3) gauge invariance of J§ through essentially the same 
argument used in Sec. IIIBT1 In fact, let U G and let U = U A denote the gauge transformation 
of U. By definition of Jq, there is a curve U(t) that solves the SD equations (j25j) . with U(0) = U 

and £7(+oo) G A/" (0) . If we define the curve U u (x; r) = U a (x;t) = A(x)U u (x; r)A(x + v) 1 , then we 

find: [7(0) = U A = U, U{+oo) G Af A = AA ), and U(r) satisfies the equation of SD by covariance 
of X7S[U\. In conclusion, although Jq cannot be expressed globally as the tensor product of SU{3) 
groups, it is nevertheless invariant under the full group of local SU{3) gauge transformation J^l. 

Note that, by definition, Jq is attached only to the middle-dimensional critical manifold A/"*- ** 
and not to the full . This is consistent with the invariance only under the SU (3) subgroup of 
SL(3, C). This also means that any section at fixed r (in a given parametrization) of the manifold 
j7o is compacJll- Therefore, gauge-fixing is not expected to be necessary to prevent numerical 
instabilities arising from gauge transformations of arbitrarily large norm. 



2. Perturbative Analysis 

We claimed that the perturbative expansion of the path integral defined on the thimble j7o 
reproduces the standard perturbation theory of QCD. In order to check this, we need to compute 
the power series in g of the integral: 

(O)o = ^[ \{dU v {x) e~ s ^O[U], Z = f \{dU v {x) e~^ u \ 

J Jo x>v J Jo XjU 

where S[U] is the action (fTU|) . 

The perturbative computation of an observable 0[A, tp, ip] at order g p involves the computation 
of integrals of the form: 

dA e S2[A]+ 9 S int [A] det (Q[ A = o]) F [ A . gi M ] q[ A = . . . . Q[A = 0; fi}- 1 ) 

J\g=0 

(28) 

In the above expression, S^fA] = | Ylx v a a(^A®(x) — /\iA^(x)) 2 , where Al is the forward lattice 
derivative. The functionals F[A;g,/j] and <Si n t[>l] are some polynomial in the field variables A. 
Note that the integral in Eq. (|28H does not include singularities in the variables A. Note also 
that the action S2 [A] does include zero modes (associated to both gauge and toronic degrees of 
freedom). These need to be regularized, as is always necessary, in lattice perturbation theory. 



11 In terms of fiber bundles, we may say that the points in J7b are sections of an S(7(3)-bundle, without being sections 
of a principal S'f7(3)-bundle. 

12 This is a vanishing (n — l)-cycle, as defined in Ref. [27t ] . 
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The expression (|28|) generates, again, the two types of terms that we have seen after Eq. (jlip . 
Those of the first type are identical to standard perturbative QCD (the argument in this case is 
even simpler, because, for g = 0, the action S2 [A] does not depend on (j,, and hence v7o(0, //) = C for 
all jj). The terms of the second type vanish for the same reason explained in relation to Eq. (I13|) . 
Hence, the perturbative series, in the expansion parameter g, of the path integrals (f2Tj) and (|2"7|) 
are identical. Note that this result is far from trivial. For example, neither the symmetries, nor 
the perturbative expansion of QCD with imaginary chemical potential [si, Ql| coincide exactly 
with those of QCD with real chemical potential. In fact, the simulations at imaginary chemical 
potential assume analyticityf^l in /i. Finally, the procedure of restricting the functional integral to 
those gaug e configurations with a positive real part of the fermionic determinant — which is known 
to fail [49] — also lacks an acceptable perturbative expansion. 



3. Relation with the Standard Approach at Zero Density 

It is interesting to check to what extent our non perturbative formulation coincides with the 
standard Wilson regularization of QCD at \x = 0. This was exactly true in the case of the model 
defined by our regularization of the action ([2]) coincide with the standard one at [i = 0, because 
the action (|2|) has one single minimum at (j) = and no further stationary points in C. 

In the case of QCD, we should first check whether the point at A = is actually a minimum 
for the QCD effective action (|19p at fi = 0. In the Appendix I A 2\ we compute the gradient and the 
Hessian of the action (|19p at A = 0. As already noticed, it is easy to check that the configuration 
A = is a stationary point. The computation of the Hessian matrix at A = requires a bit more 
work. The analytic computation is reported in Appendix IA 2\ but the final sum must be performed 
numerically. It turns out 3 that the point A = is not a minimum for periodic boundary conditions, 
but it is a minimum for (fermionic) anti-periodic boundary conditions in all directions. This is in 
agreement with the findings of [Hoi ] . In the following we will always assume this choice. 

The fact that the point A = is a local minimum, together with the observation that it is a 
global minimum in the continuum limit, justifies our approach. However, it might be interesting 
to check whether there are other (local or global) minima. For sure, there must be at least another 
local minimum, since the configuration space contains at least two disconnected components, dis- 
tinguished by the sign of the fermionic determinant and separated by a singularity of the effective 
action (|19|l 15 l. Since the minimum A = is not the only minimum, the functional integral over 
the thimble J7b attached to the point A = does not coincide with the usual functional integral 
over C, because a portion of the phase space with finite measure cannot be explored. However, 
the minimum of the component with negative fermionic determinant is most likely not a global 
minimum of the effective action, and its contribution is suppressed, as it is also confirmed by direct 



simulations 51] 



In general, the search for other minima can be done numerically by starting from random gauge 
configurations and evolving the system with the SD Eq. (125p . Since we are considering \i = 0, 
the action is real and the evolution determined by Eq. (|25p preserves the manifold C. Hence the 
evolution by SD will drive our random configuration to the local minimum of the action to which it 
is associated^!. If we repeat this procedure for a set of random configurations, we can 1) determine 
whether there are other (global or local) minima of the action (|19p besides A = 0, 2) estimate the 



Also the Taylor expansion method assumes analyticity in fi. 

This was checked for a wide range of values for the parameters g, m, r, L and T, covering those typically used in 
numerical simulations. The signs of the eigenvalues seem to depend only on the choice of the fermionic boundary 
conditions, but we could not prove this result analytically, by inspection of Eq. (|A5|) . 

Note, however, that in the complexified configuration space, the complement of the zero set of the determinant is 
connected. 

We ignore the possibility that a configuration is driven to a saddle point which is not a minimum, because this 
possibility has zero measure and hence zero probability. 
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volume of the phase space associated to the other minima and 3) compute the suppression factor 
of the non-global minima of the effective action. 

We should stress that the stationary points which are not minima (they are certainly present 
in large quantity in QCD) are irrelevant, since they represent the r — > oo limit of a set of zero 
measure in C. 

Finally, nothing that we know about non-perturbative QCD suggests that the thimble J§ might 
miss some relevant physical information. In fact, for instance, center vortices' configurations [52j 
are stationary points of the action, but not minima. Moreover, the topological sectors are not 
disconnected in lattice QCD with Wilson fermions (unless special constraints are imposed (53|). In 
particular, it is known that the iterated application of cooling transformations (which are equivalent 
to steps of SD) eventually lead any configuration with non-zero topological charge to the neutral 
topological sector (54| . Moreover, even if we consider those lattice actions that effectively separate 
the topological sectors (such as the overlap formulation [55(), it is still true that the restriction to 
the sector with zero topological charge can only introduce finite size effects to the local correlation 
functional, because no local observable can be aware of the total topological charge, in a sufficiently 
large volume. 



4- Branches of the Logarithm 

Up to now, we have used the holomorphicity of S[A] only at the critical manifold Af(°\ in order 
to deduce the properties of the Hessian matrix. The action S[A] is holomorphic in for all 
value of fx, except for a discrete set, that will be discussed in Sec. IIIICI 

Since we do not use Morse theory to justify our approach, the holomorphicity of S[A] is not 
strictly needed, besides the point A = 0. But, the insight offered by Morse theory is very important 
and we should comment on the presence of a logarithm in the effective action (1191) . The complex 
logarithmic function has two peculiarities: its imaginary part is multi-valued (or, alternatively, it 
has a cut, along one semi-axis starting from the origin of C) and it has a singularity in zero. The 
logarithm of a complex matrix has the same features for each eigenvalue of the matrix. 

The presence of the logarithm naturally leads to regard S [A] as a holomorphic function defined 
on the universal covering space X of X = SL(3,C) n \Z, where Z is the zero set of the fermionic 
determinant. Hence, the action S is holomorphic in X, which is connected and simply connected. 
This means that the thimbles attached to the stationary points of S may provide a basis for the 
homology of X , and the discussion of Sec. Ill B 31 can be repeated essentially unchanged. 

Since in numerical simulation one typically works with a parametrization of X and not of its 
covering space X, it is still necessary to check whether the multi-valuedness of the logarithmic 
function may cause some difficulties. The imaginary part of the action Si[A] is constant, by 
definition, along J§. However, it is obtained as the sum of a fermionic and a gauge part, and the 
fermionic part is itself the sum of many contributions which are not individually constant. If we 
compute the imaginary part of the logarithm in Eq. (|19|) using the prescription of the principal 
branch, we may observe jumps of 2ir. This is, however, not a problem, because the computation of 
Sj[/1] is needed only as a check of the stability of the algorithm^. Hence, any jump of a multiple 
of 2ir is acceptable. 

In principle, our set-up offers the possibility of a more elegant description of the manifold Jq 
as a true sub-manifold of the covering space X, which removes the ambiguity of 2tt completely. 

17 Note that in the case of Yang-Mills theory in 1+1 dimensions, the topological sectors affect the asymptotic space- 
time behavior of Wilson loops, also in the continuum [H(|. But in two dimensions, the Wilson loops are effectively 
global sum over the whole space-time. 

18 This is true even if we include a step of accept/reject. In this case we need to compute, besides the force 8Sr[A], 
also the value of the action However, it is always only the real part that matters. 
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In fact, any point in Jq is connected by a natural path to the manifold A/"(°), and the manifold 
AA°) is connected by a gauge transformation (which leaves all the eigenvalues of the Dirac operator 
invariant) to the point ^4 = 0. Hence, the imaginary part of the logarithm of the Dirac operator 
is well defined in any configuration in as soon as it is defined in A = 0, where it can be fixed 
conventionally. However, using this procedure to compute Sj[A] is more difficult and probably not 
justified by the wish of removing the ambiguity of 2ir. 

The singularity of the logarithm implies that the manifold Jo, defined in Sec. MI A 21 is bounded 
by the algebraic variety Z = {U : det(Q[U]) = 0}, which has complex codimension 1 (real codi- 
mension 2) in C n . This means that the curves of steepest ascent coming from AA ) at r = — oo 
may end not only at infinity (as it necessarily happens in the scalar model) but also in Z. In the 
high density regime one should expect the set Z to come very close to the point A = 0. In fact, the 
point A = actually belongs to Z for a discrete set of /i values, that become more and more dense 
in the large volume limit. This phenomenon reminds us of the fact that QCD at finite /i has really 
two — quite independent — problems: the sign problem and the problem of a high concentration of 
zero eigenvalues of the fermionic determinant near the physically interesting phase space. The two 
problems are independent, since the latter appears also where the former is absent, as in two color 



QCD 57|| . Our approach tries to address the former problem, but it is not expected to offer any 
particular advantage with respect to the latter. Nevertheless, the experience of two color QCD is 
encouraging, as it shows that much progress can be obtained in that case by tenacious algorithmic 



tunin 



C. Algorithm 

It is straightforward to adapt the algorithm described in Sec. Ill CI to the case of QCD. In this 
section we comment only on the new issues that appear in the case in QCD. In particular, the 
problem of computing the phase associated with the alignment of the tangent space of J§ with 
respect to the canonical complex volume-form is exactly the same as for the scalar theory and will 
not be discussed further. 

For convenience, we write here explicitly the main formulae. The observables that we want to 
compute can be written as: 

{0) = ±re~ iS ' [ T\dU v {x) e- s *M0[U], (29) 

and the Langevin equations are easily derived from those of SD (|25p : 

-f-U u (x;r) = -iT a (V x , u>a S\U}+Va,x,v)Uv(x-,T), (30) 

where the rj a ,x,u are random Gaussian C numbers. The procedure to project the noise vector 
into a direction tangent to Jq is exactly the same that is described in Sec. Ill CI In particular, the 
evolution equation for parallel transport of the noise vector takes the form (j, j' multi-indeces stand 
for (R/I, x, v, a)): 

±r lj (r)= Vf (r)V f V j S R , (31) 



19 Similarly to the case of two colors QCD, the drift in the Langevin equation pushes the system away from 2, which 
is important to make the problem tractable. 
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In presence of dynamical quarks, there is a further difficulty: the effective fermionic action Sf[U] 
cannot be estimated stochastically with a single extraction of pseudo-fermions, as it is usually done, 
but must be computed with sufficient precision so that the curves of SD are well defined and can 
be integrated precisely. This applies both to the evolution that makes use of Eq. (|30p and the one 
determined by Eq. (|3ip . This implies a considerable extra cost, as one Dirac inversion is necessary 
for each pseudo-fermion, and it is not clear how many pseudo-fermions will be necessary. 

An alternative could be to introduce pseudo-fermion fields and treat them similarly to the gauge 
fields. However, if we make this choice, the pseudo-fermion cannot be refreshed at every trajectory. 
Instead one should evolve them by small steps exactly like the gauge fields. It is not clear whether 
this leads to a better solution. But, it is clear that there are many possible directions in which 
one could try to improve the efficient computation of the fermionic effective action. A discussion 
of these improvements is beyond the goals of this paper. 

The procedure outlined here is certainly challenging in the case of QCD. For example, in order 
to reach the region where the quadratic approximation of the action is valid, it is necessary to apply 
as many SD iterations as necessary to suppress any nontrivial topological structure which may be 
present in the gauge configuration. The experience gained from the application of the cooling 
method 5J] suggests that the trivial sector will be reached, eventually, but it may require as many 
as 0(100) iterations. It will probably be difficult to preserve the parallelism of the noise vector r/ 
along such distance and through such nontrivial structure of the gradient flow. However, the slow 
evolution of the low modes is essentially due to the highly local nature of the smearing procedure 
which are typically used in applications that aim at preserving the low modes structure as much as 
possible. Since our goal is the opposite, we expect that techniques of Fourier acceleration \^ may 



be highly beneficial. Moreover, we expect that including a sufficient number of stout 58] or HEX 
smearing (59I ] steps in the action may reduce considerably the distance between the region where 
the simulations are performed and the region where the quadratic approximation is valid. Ideally, 
we should use an action that smears the gauge fields as much as possible, but without spoiling the 
locality of the underlying QFT. If this procedure has a chance at all, it will most probably require 
considerable tuning. 



IV. CONCLUSIONS 



We have introduced a new approach to deal with a class of sign problems that appear in lattice 
QFTs. The approach is, in principle, quite general. We have illustrated it with reference to two 
examples of QFTs with a sign problem: a scalar field theory with chemical potential and QCD at 
finite baryonic density. The former represents an ideal test bed for the method. The latter is much 
more challenging. However, it is natural to think to other possible applications, such as, e.g., the 
study of QCD with a theta term [60]. 

The approach consists in regularizing the partition function as an integral over that particular 
Lefschetz thimble which ensures a well behaved perturbative limit. There is no proof that this reg- 
ularization coincides or is physically equivalent to the standard one. Nevertheless, we have shown 
that the new formulation describes a local QFT with the correct symmetries, the correct repre- 
sentations and the correct perturbative expansion. Since the idea of universality is a fundamental 
element of our understanding of QFTs and, in particular, of QCD, this formulation does not simply 
represent a new model, but has the ambition to enable the testing of our current fundamental laws. 

Lattice QCD is too complicated to determine the exact relation between the standard formula- 
tion and the one proposed here. But, Morse theory can at least provide further evidence that the 
equivalence that we conjecture is not inconsistent. 

In this paper we have also introduced an algorithm to achieve an importance sampling of the 
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configurations in the Lefschetz thimble. It involves an elegant application of the idea of gradient 
flow to ensure that Monte Carlo updates remain in the thimble. Moreover, we have shown that 
the algorithm is protected against the obvious sources of instabilities. We certainly expect that its 
numerical application will be very challenging, especially in QCD, but we do not see any "no-go" 
obstacle that cannot be cured by a careful tuning. 

A sign problem remains, due to the relative phase between the canonical complex volume-form 
and the basis of the tangent space to the thimble. General arguments suggest that this sign 
problem should be much milder than the original one, but we have no convincing evidence of this 
yet. Moreover, computing such phase is very expensive: 0{V 2 L^) in storage and 0(V 3 ) in flops. 
Because of this, its applicability is currently restricted, at best, to very small lattices. This might 
already be very important, since the experience of the pioneering works in lattice gauge theories 
suggests that some qualitative features of the phase structure might be visible already there. 

However, the main goal of this paper is to introduce a new approach, and to prove that its 
theoretical framework is solid and that its numerical applicability is worth testing. This is necessary, 
in view of the fact that such tests (which are presently in progress) will certainly be demanding. 
On the other hand, we also hope to stimulate interest in this original approach that certainly still 
has much room for improvements. 
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Appendix A: Computation of the Hessian matrix 



In this paper we use the following notation: 

= 2sin ^ k 2 = tfi \ = sinA^ 1? = ^ V = L 3 T 

n n 

For periodic boundary conditions (which are always assumed for bosonic fields) the momenta take 
the values: 

27TT7 

Pu = -r JL n„ = 0...L„-l, (Al) 

while, in case of anti-periodic boundary conditions, the possible values are: 

P* = — ^ 2 - n„ = 0...L u -l. (A2) 
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1. The Hessian of the scalar field theory at (f> = 



Here we report the analytic computation of the Hessian matrix at (j) = 0, derived from the 
action ([5]). The Hessian matrix in configuration space reads: 



TJCLC 

x,y 



5S 



(A3) 



(2d + m 2 )5 ca S xy + 2X5 yx <j) a!X 4> CtX + A ^(((>b,x) 2 S ac 6 xy + 

b 

Sac Yl ( s v,x+i> + by,x-v) ~ 5 ac coshfi(5 y x+6 + 5 yx _ 6 ) + i sinh n £ac(S yiX+ o - * WjX _fi)» 



i>=1,3 



where e a 6 is the anti-symmetric tensor. In momentum space we have: 



Tjac r 



c.y 



e ip - x e- iq - y H^ = 0] 



<5ac m 2 + 4 sin 2 p u + 2(1 — cosh fi cos po) I — ie ac 2 sinh// sin po 

u=l,3 



(A4) 



2. The Hessian of QCD at U = 1 



Here we report the analytic computation of the Hessian matrix at A = 0, derived from the action 
j~9|) . We adopt the notation for momenta (|A1|) and (|A2|) . The covariant derivative is defined as: 



I |a=0 ' 



V x ,v,af[U] := ^/ [e iaTt ^( 



The action (1191) consists of two terms: 



S[U] = S G [U] + N f S F [U] 



S G [U] = p 
S F [U] = Ti log Q[U] 



E 



1 - -MTrU^(z) 



where the Dirac operator Q[U] is defined in (|20p . The contribution to the Hessian matrix due to 
Sq is well known: 



q( 2 ) ab („\— 1 ig(w+<5-/2) ig'(2+i > /2)v7 V7 c r/rl 



The Hessian matrix of the fermionic effective action does not seem to be available in the lit- 
erature. Hence, we report it here in some detail. We consider both periodic and anti-periodic 
boundary conditions for the fermionic fields. The following expressions are valid in both cases, but 
the substitution rule (|Alj) should be understood in the periodic case, while the rule (IA2j) holds in 
the anti-periodic case. Note that in the following the momentum variables p and k are associated to 
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fermionic degrees of freedom, while q and q' to gauge one. The derivative of the fermionic effective 
action has two contributions: 

Vw,*,aVz,u,bS F [U] = Tr {Q~ l V w ^ a V ZyV>b Q) - Tr {Q~ l {V w ^aQ)Q~ l {^ z ^ b Q)) = 

= ^ 1 Tr sc \Q X y fflw,cr,a?Jz,v,bQ)yx\ ~ ^ ^ Tr sc Q X y (V w,a,aQ)yx'Q x'y'ffl z,u,bQ) 



y'x 



The fermionic propagator reads, in momentum space: 



Qv.a ~ 3p 



m + §p 2 + if 



which is valid for both periodic and anti-periodic boundary conditions, and also at finite /i, if we 
take into account (|A1|) . (|A2|) and (fA6]). The first and second covariant derivatives of the fermion 
matrix (I20D are: 



y.z 



^ w,a,a^ z,u,bQ[U]xy — ^a,u^w,z ^ ^ T'i ) T a U u (z)e^' '° 5 X -\-u, y&x,z ^ U u (yZ)T a T b €. ^ ' ( $x—0,y&y,z^ 

In momentum space they become: 

Vz,u,bQ P , q P = 1] = ^^e^e-^V^ )6 (Q[C/ = 1])^ = 



1 



i—T b e 



iz{p-q) 



v V' 2 ' 2 



^-^A,, (L-J^T^Tae-^ + T -^T a T b e^ 



Hence, we get: 

^w,a,a^z,u,bSF[U = 1] = ^ Tr sc ^Q^ 1 ( V TO)(J , a V*,„,&<2) 

Tr sc |^(5 pfc 1 (V Wi(J>a (5)fep'(3p/fc/(V 2i ^feQ) 



+ 



fc'p 



pkp'k' 



5a,Jw,z5gb 4[(m + §j3 2 )rcosp^ - p 2 ] 

8V 2 ^ [{m + |p2)2 + p2] [( m + r £2)2 + jfe2] ' 



where 



Tip, k, a, u; m) := Tr, 



(m + -p + ^)[r(e ip<T - e ifttT ) + 7 CT (e ip,T + e~ ""HI x 
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x (m H — P 
v 2 



) + 7u(e ikv + e 



£ ^ £ l^{M p M k -pPk p )+ 

+ £ p,k £ S fa'MpMk + K,aP P kp - P a K - P v ka) + 

+ ir£ p,k £k',p(M p K + M k p u ) + ir£$£2-(Mje + M k p a ) 



and 



£ 



p.k 



,JPa 



±e 



tforr 



(m + -p 



Finally, we obtain the fermionic contribution to the Hessian in momentum space: 
1 ^^2)^(^/2)^^^ = 



o(2)oi>. 



V 

wz 

V 



E[(m + §p 2 )r cosp^ - p 2 ] 



{q<j-q v ) 



(m + §p 2 ) 2 + p 2 



T(p, fc, a, z^; m) 



[(m + §p 2 ) 2 + p 2 } [(m + §fc 2 ) 2 + k 2 ] ,, 



|fc=P+g 



(A5) 



These expressions can be generalized to the case with chemical potential p, ^ by substituting all 
the fermionic momenta {p and k, in the formulae above) as: 

Po ->■ Po + *A*. (A6) 
Eq. (|A5p is not very transparent, but can be easily computed and diagonalized numerically. 
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